package net.yegong.matrix;

import java.util.Arrays;

import org.netlib.blas.Sgemm;
import org.netlib.blas.Ssyrk;

/**
 * You must make sure these before extends this class
 * <ol>
 * <li>The matrix must store in the field <i>array</i>.
 * <li>The matrix must store in <b>column-major style</b> array. (Fortran style)
 * <li>The fields <i>nRow</i> and <i>nCol</i> must be set.
 * </ol>
 * If your matrix can't take them. Please just implement {@link SpecialMatrix}.
 * And <b>DON'T</b> extend this.
 * 
 * @author cooper
 * @see SpecialMatrix
 */
public class Matrix implements Cloneable {

	static public final MatrixType COLUMN_MAJOR = MatrixType.FORTRAN;
	static public final MatrixType ROW_MAJOR = MatrixType.C;
	static public final MatrixType FORTRAN = MatrixType.FORTRAN;
	static public final MatrixType C = MatrixType.C;

	public enum MatrixType {
		FORTRAN, C
	}

	public Matrix(float[] matrix, int nRow, int nCol, MatrixType type) {
		if (matrix.length != nRow * nCol)
			throw new MatrixDimensionException();
		this.nRow = nRow;
		this.nCol = nCol;
		this.array = matrix;
		if (type != FORTRAN) {
			array = MatrixHelper.transpose(array, nRow, nCol);
		}
	}

	public Matrix(int nRow, int nCol) {
		this.nRow = nRow;
		this.nCol = nCol;
		this.array = new float[nRow * nCol];
	}

	public Matrix(Matrix m) {
		this.nRow = m.nRow;
		this.nCol = m.nCol;
		this.array = MatrixHelper.copy(m.array);
	}

	public Matrix(SpecialMatrix m) {
		Matrix m2 = m.toMatrix();
		this.nRow = m2.nRow;
		this.nCol = m2.nCol;
		this.array = m2.array;
	}

	public final int getColumnsCount() {
		return nCol;
	}

	public final int getRowsCount() {
		return nRow;
	}

	public final int getElementsCount() {
		return nCol * nRow;
	}

	public float[] getColumnMajorArray() {
		return array;
	}

	public float[] getColumnMajorArrayCopy() {
		return MatrixHelper.copy(array);
	}

	public float[] getRowMajorArrayCopy() {
		return transpose().getColumnMajorArray();
	}

	public final float get(int i, int j) {
		return array[i + j * nRow];
	}

	public void set(int i, int j, float value) {
		array[i + j * nRow] = value;
	}

	public RowVector getRow(int index) {
		float[] array2 = new float[nCol];
		for (int i = 0; i < nCol; ++i) {
			array2[i] = array[index + i * nRow];
		}
		return new RowVector(array2);
	}

	public ColumnVector getColumn(int index) {
		float[] array2 = new float[nRow];
		System.arraycopy(array, index * nRow, array2, 0, nRow);
		return new ColumnVector(array2);
	}

	public void setRow(int index, RowVector row) {
		if (nCol != row.nCol)
			throw new MatrixDimensionException();
		for (int i = 0; i < nCol; ++i) {
			array[index + i * nRow] = row.array[i];
		}
	}

	public void setRow(int index, float value) {
		for (int i = 0; i < nCol; ++i) {
			array[index + i * nRow] = value;
		}
	}

	public void setColumn(int index, ColumnVector col) {
		if (nRow != col.nRow)
			throw new MatrixDimensionException();
		System.arraycopy(col.array, 0, array, index * nRow, nRow);
	}

	public void setColumn(int index, float value) {
		Arrays.fill(array, index * nRow, (index + 1) * nRow, value);
	}

	public float[] getRowArray(int index) {
		float[] array2 = new float[nCol];
		for (int i = 0; i < nCol; ++i) {
			array2[i] = array[index + i * nRow];
		}
		return array2;
	}

	public float[] getColumnArray(int index) {
		float[] array2 = new float[nRow];
		System.arraycopy(array, index * nRow, array2, 0, nRow);
		return array2;
	}

	public Matrix transpose() {
		float[] array2 = MatrixHelper.transpose(array, nCol, nRow);
		return new Matrix(array2, nCol, nRow, FORTRAN);
	}

	public void transposeMe() {
		int tmp = nCol;
		nCol = nRow;
		nRow = tmp;
		array = MatrixHelper.transpose(array, nRow, nCol);
	}

	public Matrix add(Matrix m) {
		if (nRow != m.nRow || nCol != m.nCol)
			throw new MatrixDimensionException();
		Matrix ret = clone();
		int len = array.length;
		for (int i = 0; i < len; ++i) {
			ret.array[i] += m.array[i];
		}
		return ret;
	}

	public Matrix minus(Matrix m) {
		if (nRow != m.nRow || nCol != m.nCol)
			throw new MatrixDimensionException();
		Matrix ret = clone();
		int len = array.length;
		for (int i = 0; i < len; ++i) {
			ret.array[i] -= m.array[i];
		}
		return ret;
	}

	public Matrix add(Matrix m, float alpha, float beta) {
		if (nRow != m.nRow || nCol != m.nCol)
			throw new MatrixDimensionException();
		Matrix ret = new Matrix(nRow, nCol);
		int len = array.length;
		for (int i = 0; i < len; ++i) {
			ret.array[i] = alpha * array[i] + beta * m.array[i];
		}
		return ret;
	}

	public Matrix multiply(float alpha) {
		Matrix ret = clone();
		int len = array.length;
		for (int i = 0; i < len; ++i) {
			ret.array[i] *= alpha;
		}
		return ret;
	}

	public void scale(float alpha) {
		int len = array.length;
		for (int i = 0; i < len; ++i) {
			array[i] *= alpha;
		}
	}

	public SymmetricMatrix leftMultiplyMyTranspose() {
		float[] array2 = new float[nRow * nRow];
		Ssyrk.ssyrk("U", "N", nRow, nCol, 1, array, 0, nRow, 0, array2, 0, nRow);
		for (int i = 0; i < nRow; ++i) {
			for (int j = i + 1; j < nRow; ++j) {
				array2[j + i * nRow] = array2[i + j * nRow];
			}
		}
		SymmetricMatrix sym = new SymmetricMatrix(array2, nRow, FORTRAN);
		return sym;
	}

	public SymmetricMatrix rightMultiplyMyTranspose() {
		float[] array2 = new float[nCol * nCol];
		Ssyrk.ssyrk("U", "T", nCol, nRow, 1, array, 0, nRow, 0, array2, 0, nCol);
		for (int i = 0; i < nCol; ++i) {
			for (int j = i + 1; j < nCol; ++j) {
				array2[j + i * nCol] = array2[i + j * nCol];
			}
		}
		SymmetricMatrix sym = new SymmetricMatrix(array2, nCol, FORTRAN);
		return sym;
	}

	public float trace() {
		return getDiag().trace();
	}

	/**
	 * @param m Matrix M
	 * @return A * M
	 */
	public Matrix leftMultiply(Matrix m) {
		if (m instanceof ColumnVector)
			return m.rightMultiply(this);
		if (nCol != m.nRow)
			throw new MatrixDimensionException();
		float[] array2 = new float[nRow * m.nCol];
		Sgemm.sgemm("N", "N", nRow, m.nCol, nCol, 1, array, 0, nRow, m.array, 0, m.nRow, 0, array2, 0, nRow);
		return new Matrix(array2, nRow, m.nCol, FORTRAN);
	}

	/**
	 * @param m Matrix M
	 * @return M * A
	 */
	public Matrix rightMultiply(Matrix m) {
		return m.leftMultiply(this);
	}

	public void reserveUpperTriangle() {
		int k = Math.min(nRow, nCol);
		for (int i = 0; i < k; ++i) {
			Arrays.fill(array, i + 1 + i * nRow, (i + 1) * nRow, 0);
		}
	}

	public Matrix getUpperTriangle() {
		Matrix ret = clone();
		ret.reserveUpperTriangle();
		return ret;
	}

	public void reserveLowerTriangle() {
		for (int i = 0; i < nCol; ++i) {
			Arrays.fill(array, i * nRow, Math.min(i, nRow) + i * nRow, 0);
		}
	}

	public Matrix getLowerTriangle() {
		Matrix ret = clone();
		ret.reserveLowerTriangle();
		return ret;
	}

	public DiagonalMatrix getDiag() {
		int len = Math.min(nRow, nCol);
		float[] diag = new float[len];
		for (int i = 0; i < len; ++i) {
			diag[i] = array[i + i * nRow];
		}
		return new DiagonalMatrix(diag, nRow, nCol);
	}

	/**
	 * This method will solve the equation A*x = b
	 * 
	 * @param b Vector b
	 * @return Vector x.
	 */
	public ColumnVector solve(ColumnVector b) {
		return MatrixHelper.solve(this, b);
	}

	/**
	 * This method will solve the equation x*A = b
	 * 
	 * @param b Vector b
	 * @return Vector x.
	 */
	public RowVector solve(RowVector b) {
		return MatrixHelper.solve(this, b);
	}

	public Matrix toMatrix() {
		return this;
	}

	public float max() {
		float max = Float.NEGATIVE_INFINITY;
		for (float d : array)
			if (d > max)
				max = d;
		return max;
	}

	public ColumnVector maxOfRow() {
		float[] max = new float[nRow];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			float aMax = Float.NEGATIVE_INFINITY;
			for (int j = 0; j < nCol; ++j) {
				tmp = array[i + j * nRow];
				if (tmp > aMax)
					aMax = tmp;
			}
			max[i] = aMax;
		}
		return new ColumnVector(max);
	}

	public RowVector maxOfColumn() {
		float[] max = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			float aMax = Float.NEGATIVE_INFINITY;
			for (int j = 0; j < nRow; ++j) {
				tmp = array[j + i * nRow];
				if (tmp > aMax)
					aMax = tmp;
			}
			max[i] = aMax;
		}
		return new RowVector(max);
	}

	public float maxAbsolute() {
		float max = Float.NEGATIVE_INFINITY;
		float abs;
		for (float d : array) {
			abs = Math.abs(d);
			if (abs > max)
				max = abs;
		}
		return max;
	}

	public ColumnVector maxAbsoluteOfRow() {
		float[] max = new float[nRow];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			float aMax = Float.NEGATIVE_INFINITY;
			for (int j = 0; j < nCol; ++j) {
				tmp = Math.abs(array[i + j * nRow]);
				if (tmp > aMax)
					aMax = tmp;
			}
			max[i] = aMax;
		}
		return new ColumnVector(max);
	}

	public RowVector maxAbsoluteOfColumn() {
		float[] max = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			float aMax = Float.NEGATIVE_INFINITY;
			for (int j = 0; j < nRow; ++j) {
				tmp = Math.abs(array[j + i * nRow]);
				if (tmp > aMax)
					aMax = tmp;
			}
			max[i] = aMax;
		}
		return new RowVector(max);
	}

	public float min() {
		float min = Float.POSITIVE_INFINITY;
		for (float d : array)
			if (d < min)
				min = d;
		return min;
	}

	public ColumnVector minOfRow() {
		float[] min = new float[nRow];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			float aMin = Float.POSITIVE_INFINITY;
			for (int j = 0; j < nCol; ++j) {
				tmp = array[i + j * nRow];
				if (tmp < aMin)
					aMin = tmp;
			}
			min[i] = aMin;
		}
		return new ColumnVector(min);
	}

	public RowVector minOfColumn() {
		float[] min = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			float aMin = Float.POSITIVE_INFINITY;
			for (int j = 0; j < nRow; ++j) {
				tmp = array[j + i * nRow];
				if (tmp < aMin)
					aMin = tmp;
			}
			min[i] = aMin;
		}
		return new RowVector(min);
	}

	public float minAbsolute() {
		float min = Float.POSITIVE_INFINITY;
		float abs;
		for (float d : array) {
			abs = Math.abs(d);
			if (abs < min)
				min = abs;
		}
		return min;
	}

	public ColumnVector minAbsoluteOfRow() {
		float[] min = new float[nRow];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			float aMin = Float.POSITIVE_INFINITY;
			for (int j = 0; j < nCol; ++j) {
				tmp = Math.abs(array[i + j * nRow]);
				if (tmp < aMin)
					aMin = tmp;
			}
			min[i] = aMin;
		}
		return new ColumnVector(min);
	}

	public RowVector minAbsoluteOfColumn() {
		float[] min = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			float aMin = Float.POSITIVE_INFINITY;
			for (int j = 0; j < nRow; ++j) {
				tmp = Math.abs(array[j + i * nRow]);
				if (tmp < aMin)
					aMin = tmp;
			}
			min[i] = aMin;
		}
		return new RowVector(min);
	}

	public float sum() {
		float sum = 0;
		for (float d : array) {
			sum += d;
		}
		return sum;
	}

	public ColumnVector sumOfRow() {
		float[] sum = new float[nRow];
		for (int i = 0; i < nRow; ++i) {
			for (int j = 0; j < nCol; ++j) {
				sum[i] += array[i + j * nRow];
			}
		}
		return new ColumnVector(sum);
	}

	public RowVector sumOfColumn() {
		float[] sum = new float[nCol];
		for (int i = 0; i < nCol; ++i) {
			for (int j = 0; j < nRow; ++j) {
				sum[i] += array[j + i * nRow];
			}
		}
		return new RowVector(sum);
	}

	public float mean() {
		return sum() / array.length;
	}

	public ColumnVector meanOfRow() {
		float[] mean = new float[nRow];
		for (int i = 0; i < nRow; ++i) {
			for (int j = 0; j < nCol; ++j) {
				mean[i] += array[i + j * nRow];
			}
		}
		for (int i = 0; i < nRow; ++i) {
			mean[i] /= nCol;
		}
		return new ColumnVector(mean);
	}

	public RowVector meanOfColumn() {
		float[] mean = new float[nCol];
		for (int i = 0; i < nCol; ++i) {
			for (int j = 0; j < nRow; ++j) {
				mean[i] += array[j + i * nRow];
			}
		}
		for (int i = 0; i < nCol; ++i) {
			mean[i] /= nRow;
		}
		return new RowVector(mean);
	}

	/**
	 * @return population standard deviation
	 */
	public float stdev(float mean) {
		float stdev = 0;
		for (float d : array) {
			stdev += (d - mean) * (d - mean);
		}
		return (float) Math.sqrt(stdev / array.length);
	}

	/**
	 * @return population standard deviation
	 */
	public ColumnVector stdevOfRow(ColumnVector means) {
		float[] mean = means.array;
		float[] stdev = new float[nRow];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			for (int j = 0; j < nCol; ++j) {
				tmp = array[i + j * nRow];
				stdev[i] += (tmp - mean[i]) * (tmp - mean[i]);
			}
			stdev[i] = (float) Math.sqrt(stdev[i] / nCol);
		}
		return new ColumnVector(stdev);
	}

	/**
	 * @return population standard deviation
	 */
	public RowVector stdevOfColumn(RowVector means) {
		float[] mean = means.array;
		float[] stdev = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			for (int j = 0; j < nRow; ++j) {
				tmp = array[j + i * nRow];
				stdev[i] += (tmp - mean[i]) * (tmp - mean[i]);
			}
			stdev[i] = (float) Math.sqrt(stdev[i] / nRow);
		}
		return new RowVector(stdev);
	}

	public float sumSquare() {
		float sum = 0.0f;
		for (float d : array) {
			sum += d * d;
		}
		return sum;
	}

	public ColumnVector sumSquareOfRow() {
		float[] sum2 = new float[nCol];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			for (int j = 0; j < nCol; ++j) {
				tmp = array[i + j * nRow];
				sum2[i] += tmp * tmp;
			}
		}
		return new ColumnVector(sum2);
	}

	public RowVector sumSquareOfColumn() {
		float[] sum2 = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			for (int j = 0; j < nRow; ++j) {
				tmp = array[j + i * nRow];
				sum2[i] += tmp * tmp;
			}
		}
		return new RowVector(sum2);
	}

	/**
	 * @return sample standard deviation
	 */
	public float stdev() {
		float mean = mean();
		float stdev = 0;
		for (float d : array) {
			stdev += (d - mean) * (d - mean);
		}
		return (float) Math.sqrt(stdev / (array.length - 1));
	}

	/**
	 * @return sample standard deviation
	 */
	public ColumnVector stdevOfRow() {
		float[] mean = meanOfRow().array;
		float[] stdev = new float[nRow];
		float tmp;
		for (int i = 0; i < nRow; ++i) {
			for (int j = 0; j < nCol; ++j) {
				tmp = array[i + j * nRow];
				stdev[i] += (tmp - mean[i]) * (tmp - mean[i]);
			}
			stdev[i] = (float) Math.sqrt(stdev[i] / (nCol - 1));
		}
		return new ColumnVector(stdev);
	}

	/**
	 * @return sample standard deviation
	 */
	public RowVector stdevOfColumn() {
		float[] mean = meanOfColumn().array;
		float[] stdev = new float[nCol];
		float tmp;
		for (int i = 0; i < nCol; ++i) {
			for (int j = 0; j < nRow; ++j) {
				tmp = array[j + i * nRow];
				stdev[i] += (tmp - mean[i]) * (tmp - mean[i]);
			}
			stdev[i] = (float) Math.sqrt(stdev[i] / (nRow - 1));
		}
		return new RowVector(stdev);
	}

	public float det() {
		LUDecomposition lu = new LUDecomposition(this);
		return lu.det();
	}

	public int rank() {
		QRDecomposition qr = new QRDecomposition(this);
		return qr.rank();
	}

	/**
	 * Please notice that this method compute <b>2-norm</b> of the matrix. Other
	 * norm is not included in this version. If you think it's necessary, let me
	 * know please.
	 * 
	 * @return The 2-norm of the matrix
	 */
	public float norm() {
		SingularValueDecomposition svd = new SingularValueDecomposition(this);
		return svd.norm();
	}

	/**
	 * Please notice that this method compute cond of the matrix using
	 * <b>2-norm</b>. Other norm is not included in this version. If you think
	 * it's necessary, let me know please.
	 * 
	 * @return The cond of the matrix
	 */
	public float cond() {
		SingularValueDecomposition svd = new SingularValueDecomposition(this);
		return svd.cond();
	}

	public Matrix rowSlice(int start, int length) {
		float[] array2 = new float[nCol * length];
		for (int j = 0; j < nCol; ++j) {
			System.arraycopy(array, start + j * nRow, array2, j * length, length);
		}
		if (length == nCol)
			return new SquareMatrix(array2, length, FORTRAN);
		return new Matrix(array2, length, nCol, FORTRAN);
	}

	public Matrix columnSlice(int start, int length) {
		float[] array2 = new float[nRow * length];
		System.arraycopy(array, start * nRow, array2, 0, nRow * length);
		if (length == nRow)
			return new SquareMatrix(array2, length, FORTRAN);
		return new Matrix(array2, nRow, length, FORTRAN);
	}

	public SymmetricMatrix toSymmetric() {
		return new SymmetricMatrix(array, nRow, FORTRAN);
	}

	public SquareMatrix toSquare() {
		return new SquareMatrix(array, nRow, FORTRAN);
	}

	public RowVector toRowVector() {
		return new RowVector(array);
	}

	public ColumnVector toColumnVector() {
		return new ColumnVector(array);
	}

	@Override
	public Matrix clone() {
		float[] array2 = new float[nRow * nCol];
		System.arraycopy(array, 0, array2, 0, nRow * nCol);
		return new Matrix(array2, nRow, nCol, FORTRAN);
	}

	@Override
	public String toString() {
		StringBuilder buf = new StringBuilder();
		for (int i = 0; i < getRowsCount(); ++i) {
			for (int j = 0; j < getColumnsCount(); ++j) {
				buf.append(array[i + j * nRow]);
				buf.append(", ");
			}
			buf.append("\n");
		}
		return buf.toString();
	}

	int nRow, nCol;
	protected float[] array;
}
